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Abstract 


In this paper, we present numerical methods to implement the probabilistic repre¬ 
sentation of third kind (Robin) boundary problem for the Laplace equations. The 
solution is based on a Feynman-Kac formula for the Robin problem which employs 
the standard reflecting Brownian motion (SRBM) and its boundary local time aris¬ 
ing from the Skorohod problem. By simulating SRBM paths through Brownian 
motion using Walk on Spheres (WOS) method, approximation of the boundary lo¬ 
cal time is obtained and the Feynman-Kac formula is calculated by evaluating the 
average of all path integrals over the boundary under a measure defined through 
the local time. Numerical results demonstrate the accuracy and efficiency of the 
proposed method for finding a local solution of the Laplace equations with Robin 
boundary conditions. 


Key words: Skorohod problem, boundary local time, Feynman-Kac formula. 
Reflecting Brownian Motion, Brownian motion, Laplace equation, WOS, Robin 
boundary problem 


1 Introduction 


Partial differential equations (PDEs) have been widely used to describe a 
variety of phenomena such as electrostatics, electrodynamics, fluid flow or 
quantum mechanics. Traditionally, hnite difference, hnite element and bound¬ 
ary element methods are the mainstream numerical approaches to solve the 
PDEs. Recently, using the Feynman-Kac formula |5] |6] [7] which connects solu¬ 
tions of differential equations of diffusion and heat flow and random processes 
of Brownian motions, numerical methods based on random walks or Monte 
Carlo diffusions have been explored for solving parabolic and elliptic PDEs 

niza. 

The Feynman-Kac formula represents the solutions of parabolic and ellip¬ 
tic PDEs as the expectation functionals of stochastic processes (specihcally 
Brownian motions), and conversely, the probabilistic properties of diffusion 
processes can be obtained through investigating related PDEs characterized 
by corresponding generators [22]. The formula involves the path integrals of the 
diffusion process starting from an arbitrarily prescribed location, and this en¬ 
ables us to hnd a local numerical solution without constructing space and time 
meshes as in traditional deterministic numerical methods mentioned above, 
which incur expensive costs in high dimensions. In many applications it is 
also of practical importance and necessity to seek a local solution of PDEs at 
some interested points. If the sample paths of a diffusion process are simulated, 
then by computing the average of path integrals we can obtain approxima- 




tions to the exact solutions of the PDEs. For second order elliptic PDFs with 
Dirichlet and Neumann boundaries, the average of path integrals is reduced 
to the average of boundary integrals under certain measure where the detailed 
trajectories of the diffusion process have no effect on the averages except the 
hitting locations on the boundaries. 


Simulations of diffusion paths can be done by random walks methods I2IISI 
n IE] either on lattice or in continuum space. In some cases such as for the 
Poisson equation, the Feynman-Kac formula has a pathwise integral requiring 
the detailed trajectory of each path. Moreover, one may need to adopt ran¬ 
dom walks on a discrete lattice in order to incorporate inhomogeneous source 
terms. As for the continuum space approach, the Walk on Spheres (WOS) 
method is preferred where the path of diffusion process within the domain 
does not appear in the Feynman-Kac formula. For both approaches, the ge¬ 
ometry of the boundaries need special care for accurate results [H]. In our 
previous work on Laplace equation with Neumann boundary conditions [4], 
we proposed a numerical method to simulate the standard reflecting Brown¬ 
ian motion (SRBM) path using WOS and obtained the boundary local time 
of the SRBM. As a result, a local numerical solution of the PDE is achieved 
by using the Feynman-Kac formula. Other literatures |9] |T0] |T3] [H] have also 
explored similar problems. Especially, in [13] schemes based on the WOS, Eu¬ 
ler schemes and kinetic approximations are proposed to treat inhomogeneous 
Neumann problems. It turns out that the pointwise resolution is much harder 
due to the choice of the truncation of time. However, the local time was not 
handled explicitly in [13]. On the other hand, Monte Carlo simulations were 
discussed in [T3] where the positive part of the boundary needs to be identified 
hrst. In this paper, following [1] we continue the use of SRBM to solve Robin 
boundary problems for the Laplace operator, which has many applications in 
heat transfer and impedance tomograph. Our goal again is to obtain a local 
approximation to the exact solution of the Robin problem. 


The rest of paper is organized as follows. Firstly, the Skorohod problem is in¬ 
troduced in section 2, where both the concepts of standard reflecting Brownian 
motion and boundary local time will be reviewed briefly. This lays the founda¬ 
tion for the underlying diffusion process of the Robin boundary problem and 
the sampling of the diffusion paths. Secondly, an overview of the Feynamn- 
Kac formula is given in section 3. Thirdly, the probabilistic representation of 
the solution for the Robin boundary value problem proposed in mm is dis¬ 
cussed in section 4, and we will see the relation between the Neumann and 
Robin problems and gain a new perspective. Section 5 presents the numerical 
approaches and test results. Finally, conclusions and future work are given in 
section 6. 


2 Skorohod problem, SRBM and boundary local time 


Assume that D is a domain with a boundary in R^. The generalized Sko¬ 
rohod problem is stated as follows: 


Definition 1 Let f G ^([O, cxo), a continuous function from [0, oo] to R?. 
A pair {^t,Lt) is a solution to the Skorohod equation S{f-,D) if 

(1) f is continuous in D; 

(2) L(t) is a nondecreasing function which increases only when f G dD, 
namely, 

L{t)= f lanimWs); (2.1) 

Jo 

(3) The Skorohod equation holds: 

S(f;D): at) = m - (2.2) 

where n{x) denotes the outward unit normal vector at x E dD. 


The Skorohod problem was hrst studied in [T] by A.V. Skorohod in addressing 
the construction of paths for diffusion processes with boundaries, which results 
from the instantaneous reflection behavior of the processes at the boundaries. 
Skorohod presented the result in one dimension in the form of an Ito integral 
and Hsu [12] later extended the concept to d-dimensions (d > 2). 


In the simple case that D = [0,cxd), the solution to the Skorohod problem 
uniquely exists and can be explicitly given by 




fit), if t<T-, 

fit)- inf /(s), ift>T; 


(2.3) 


where r = inf {f > 0 : /(f) < 0}. In general, solvability of the Skorohod prob¬ 
lem is closely related to the smoothness of the domain D. For higher dimen¬ 
sions, the existence of fl2.2p is guranteed for domains while uniqueness can 
be acheived for a domain by assuming the convexity for the domain [T5] . 
Later, it was shown by Lions and Sznitman [T6| that the constraints on D can 
be relaxed to some locally convex properties. 


Next we introduce the concept of SRBM and boundary local time which 
play important roles in solving Robin boundary problem by probabilistic ap¬ 
proaches. 

Suppose that fit) is a standard Brownian motion (SBM) starting at x E D 
and iXt, Lt) is the solution to the Skorohod problem S'(/; D), then Xt will be 


the standard reflecting Brownian motion (SRBM) on D starting at x. Becanse 
the transition probability density of the SRBM satisfies the same parabolic 
differential eqnation as that by a BM, a sample path of the SRBM can be 
simnlated simply as that of the BM within the domain. However, the zero 
Nenmann bonndary condition for the density of SRBM implies that the path 
be pnshed back at the bonndary along the inward normal direction whenever 
it attempts to cross the latter. The fnll constrnction of a SRBM from a SBM 
can be fonnd in onr previons work [1]. 

The bonndary local time is not an independent process bnt associated with 
SRBM Xt and defined by 


Lit) = lim 
^ ' 6-|>0 


t,lDXXs)ds 

6 


(2.4) 


where is a strip region of width e containing dD and G D. Here Lt is 
called the local time of X^, a notion invented by P. Levy [23]. This limit exists 
both in and P^-a.s. for any x E D. 

It is obvions that Lt measures the amount of time that the standard reflect¬ 
ing Brownian motion Xt spends in a vanishing neighborhood of the boundary 
within the time period [0, t]. Besides, it is the unique continuous nondecreasing 
process that appears in the Skorohod equation. An interesting part of (12.4p is 
that the set {t G R+ : Xt G dD} has a zero Lebesgue measure while the so¬ 
journ time of the set is nontrivial [23]. This concept is not just a mathematical 
one but also has physical relevance in understanding the “crossover exponent” 
associated with “renewal rate” in modern renewal theory im. 

In [12], an alternative explicit form of the local time was found, 

= tldD{Xs)VTs, (2.5) 

V 2 JO 

where the the right-hand side of fl2.5p is understood as the limit of 


n—1 


i=l 

where A = {Aj} is a partition of the interval [0, t] and each Aj is an element in 
A. fl2.4p and fl2.5l) provide us different ways to approximate local time and in 
[1], it was found that fl2.4p yields better approximations in Neumann problem 
than (12.5p . Therefore, in this paper, we will also choose (12.4p as the approach 
to estimate the local time here. 



maxIgoiXs 

sG Aj 


Aj|, max|Ai| —)■ 0, 


( 2 . 6 ) 













3 A Feynman-Kac formula 


The Feynman-Kac formula named after Richard Feynman and Mark Kac, 
establishes a link between PDFs and stochastic processes. It hrst arose in the 
potential theory for Schodinger equations, leading to a profound reformulation 
of the quantum mechanics by the means of path integrals. Later, the formula 
also hnds its applications in mathematical finance, where the probabilistic and 
the PDF representations in derivative pricing are connected. 


Let us first look at the Dirichlet problems. Given a domain D <Z with a 
boundary dD, 


Lu{x) — c{x)u{x) = f{x), X E D 

u{x) = X G dD ’ 


(3.1) 


where the operator L = -- ^2 


*j=i 


32 d d 

7 -— —r — bAx)—— and both the 

dx^dx^ ^ ^ ^dx^ 


coefficients in L and c{x) are Lipschitz continuous and bounded. 


The Feynman-Kac formula in this case [18] represents the solution to fl3.ip in 
terms of an Ito diffusion process XA^o), 


TD 


r^D 


u{x) = E^[J^ f{Xt)exp^J^ c{Xs)dsjdt] +E'"[(j){Xrjy)exp^J^ c(Xs)(is|], 

(3.2) 

with td = inf{t : Xt G dD} and XAu) is defined by 


dXt = b(Xt)dt + a{Xt)dBt, 


(3.3) 


where Bt is the Brownian motion and [a^] 


1 

2 


a{x)a^{x), [bij] 


= b. 


The expectation E^ is an integration with respect to a measure Px taken over 
all sample paths Xt=o{Lj) = x, thus fl3.2p is a representation of a solution of 
Dirichlet problem in the form of functional integral. Moreover, fl3.2p is obtained 
by killing process Xt at a stopping time td at which Xt will be absorbed on the 
boundary. If c{x) > 0, then the function c{x) can be interpreted as the killing 
rate ■ R should be pointed out that fl3.2p is equivalent to the formulation of 
weak solution and it is a classical solution as well if some smoothness conditions 
are satisfied. 


The Feynman-Kac formula above offers a method for solving certain PDFs 
by simulating random paths of a stochastic process. Conversely, an important 
class of expectations of random processes can be computed by deterministic 
methods. For the Neumann boundary condition, a similar formula was derived 
by Hsu 113 for the Poisson equation, which is in the form of a functional 
integral based on the boundary local time introduced in section 2. In this 







case, though the Feynman-Kac formula remains in a similar form, it should 
be understood as a path integral over the stochastic process Lt associated with 
the standard reflecting Brownian motion. 


4 Robin boundary value problem 


We focus on Robin boundary value problem for the time-independent Schrodinger 
equation. 


1 


Au + qu = 0, in D; 

^ — cu = f, on dD. 
on 


(4.1) 


A generalization of the Feynman-Kac formula of section 3 in [2] gives a prob- 
ablistic solution of fl4.ip as follows. 


u(x) = 


eg(t)ecit)f(Xt)dLt \ , 


(4.2) 


where Xt is a SRBM starting at x. The term Feynman-Kac functional eq{t), 
also appeared in the Neumann problem [T2], is defined as 


(4.3) 


eq{t)=exp [ q{Xs)ds 
Uo 


and a second functional is introduced for the Robin boundary problem, for 
C e MdD) 


Ccit) = exp 


Uo 


ci^Xg^dLs 


Using these two functionals, we have, 

ft 


uix) = 


/o 


exp 


/o 


{q{Xs)ds + c{Xs)dL,) 


f{Xt)dL, . 


(4.4) 


(4.5) 


L{t) ~ / lDXXs)ds, 

e Jo 


Recalling the dehnition of the local time in fl2.4l) . we have the following ap¬ 
proximation 

(4.6) 

(4.7) 


thus. 


dL{s) ^ -lDSXs)ds. 


Therefore, fl4.5p can be modified as 

rt / 1 


u{x) ^ E^ exp (g(X,) + ^c(X,)/z,,(X,)) ds] /(W)dLi|, (4.8) 















It can also be shown that as e goes to zero, fl4.8p converges to fl4.5p nniformly 
on D. 

As (14.Sp resembles the Feyman-Kac formnla for the Nenmann problem with a 
modihed q{x) |1], it indicates a connection between the Robin and the Nen¬ 
mann problems, namely, we may introdnce 

q^{x) = q{x) + ^c{x)Id,{x), (4.9) 

then, the Robin bonndary problem fl4.ip can be viewed as a limiting case 
(e —)■ 0) of Neumann problems 

{ ^An -h q^u 

du 
dn 


0, in D; 

/, on dD. 


(4.10) 


5 Numerical approach and results 


In the present work, we only consider the case of the Laplace equation where 
g = 0 in (I4.10p . From (14.bp . 

u{x) = F;" , (5.1) 

where represents the standard reflecting Brownian motion. For the sake of 
computer simulation, the time period is truncated into [0,T] to produce an 
approximation for u{x), i.e.. 


u{x) = 




(5.2) 


Next we will give a general description on the realization of SRBM paths and 
the calculation of the corresponding local time, as implemented in |1]. A SRBM 
path can be constructed by pulling back a BM path back onto the boundary 
whenever it runs out of the domain. Specihcally, a SRBM path behaves exactly 
the same way as a BM which is simulated by the WOS method. 


5.1 Simulating SRBM by the method of Walk on Spheres (WOS) 
• Method of WOS for Brownian paths 










Random walk on spheres (WOS) method was hrst proposed by Miiller |^, 
which can solve the Dirichlet problem for the Laplace operator efficiently 

HP! . 


To illustrate the WOS method for the Dirichlet problem (13.Ih . let us consider 
the Laplace equation again where / = 0, Qij = Sij and 6^ = 0 in fl3.ip and the 
ltd diffusion is then simply the standard Brownian motion with no drift. The 
solution to the Laplace equation can be rewritten in terms of a measure 
dehned on the boundary dD, 

u{x) = (5.3) 

JdD 

where is the harmonic measure dehned by 

^l{F) = P^{Xr^eF},F ^dD,xeD. (5.4) 


It can be shown easily that the harmonic measure is related to the Green’s 
function g{y,x) for the domain with a homogeneous boundary condition |20j . 


i.e., 


as follows 


^ 9 {x,y) = 6{x 

-y), xeD, 

(5.5) 

o' 

X e dD ’ 

p(x,y) = 

dg{x,y) 

dn,. 

(5.6) 


If the starting point x of a Brownian motion is at the center of a ball, the 
probability of the BM exiting a portion of the boundary of the ball will be 
proportional to the portion’s area. Therefore, sampling a Brownian path by 
drawing balls within the domain can signihcantly reduce the path sampling 
time. To be specihc, given a starting point x inside the domain D, we simply 
draw a ball of largest possible radius fully contained in D and then the next 
location of the Brownian path on the surface of the ball can be sampled, 
using a uniform distribution on the sphere, say at xi. Treat xi as the new 
starting point, draw a second ball fully contained in D, make a jump from xi 
to X 2 on the surface of the second ball as before. Repeat this procedure until 
the path hits a absorption e-shell of the domain (see Fig. 2) |5]. When this 
happens, we assume that the path has hit the boundary dD (see Fig. 1(a) for 
an illustration). 

Now we can dehne an estimator of fl3.2p with c = 0 by 

1 ^ 

~ ( 5 - 7 ) 

i=i 

where N is the number of Brownian paths sampled and x* is the hrst hitting 
point of each path on the boundary. To speed up the WOS process, maximum 





(a) WOS within the domain 


(b) WOS (with a maximal 
step size for each jump) within 
the domain 


Fig. 1. Walk on Spheres method 

possible size of the sphere for each step would allow faster hrst hitting on the 
boundary. 


• WOS and RBM 


For the reflecting boundary, we will construct a strip region around the bound¬ 
ary (see Fig. 2) and allow the process Xt to move according to the law of BM 
continuously. Before the path enters the strip region, the radius of WOS is 
chosen to be of a maximum possible size less than the distance to the bound¬ 
ary. Once the particle is in the strip region, the radius of the WOS sphere is 
hxed at a constant Xx (or 2Xx, see Fig. 3). With this approach, according to 
the dehnition fl2.4p . the local time may be interpreted as 


dL{t) 


ILIdAxMs 


(6.8) 


which is 


dL{t) 


IL lD,(.X.)ds 


= (ni, - Ki,-.) 


{AxY 


(5.9) 


e " " * 3e 

given a prehxed constant Ax in the strip region and ntj be the cumulative 
steps that path stays within the e-region from the begining until time tj (see 
Remark below for dehnition). Notice that only those steps where the path of 
Xt remains in the e-region will contribute to rit^ because the SRBM may he 
out of the e-region at other steps. More details can be found in [1], where 
the same construction is applied for the Neumann boundary value problem. 
One may refer to Fig. 3 for an illustration of the behavior of path near the 
boundary. 


Remark 2 Occupation time of SRBM Xt in the numerator of / (5.8I) was cal¬ 
culated in terms of that of BM sampled by the walks on spheres. Notice here 
that within the e-region, the radius of the WOS may he Ax or 2Ax, which 
implies that the corresponding elapsed time of one step for local time could be 
(Aa;)^/3 or (2Aa:)^/3. The latter is four times bigger than the former. But if 
we absorb the factor 4 into nt, (15.91) still holds. In practical implementation, 












£‘Shell 



Fig. 2. A e-region for a bounded domain in 



Fig. 3. WOS in the e-region. At point xi, BM path hrst hits the e-region. By WOS 
with a prefixed radius Ax, the path continues moving subsequently to X 2 where the 
distance to the boundary is less than Ax. Enlarge the radius to 2Ax, the path then 
have a probability to run out of the domain to X 3 . Pull back to the closest point X 4 
on the boundary, record </>(x4) and continue WOS-sampling starting at X4. 

we treat rit as a vector of entries of increasing value, the increment of each 
component of Ut over the previous one after each step of WOS will be 0, 1 or 
4, corresponding to the scenarios that Xt is out of the e-region, in the e-region 
while sampled on the sphere of a radius Ax, or in the e-region while sampled 
on the sphere of a radius 2Ax, respectively. 

Robin boundaries represent a general form of an insulating boundary con¬ 
dition for convection-diffusion equations where c{x) stands for the positive 
diffusive coefficients. For our numerical test, we will consider two cases: a 
positive constant c and a positive function c{x). 


5.2 Numerical Tests 


The numerical approximations obtained are compared to the true solutions on 
a selected circle and a line segment, respectively, for the following three test 
domains in R^: 




(1) A cube centered at the origin with a length 2; 

(2) A sphere centered at the origin with a radius 1; 

(3) An ellipsoid centered at the origin with axial lengths [3, 2, 1], 


The location of the circle is given by 

{{x.y^zY = (r cos 6^1 sin 02 , r sin 01 sin 6 ^ 2 , r cos 6^2 )^} (5.10) 

with r = 0.6, 0i = 0 : fc ■ 27r/30 : 27r, 02 = vr/d with fc = 1,..., 15. While the 
line segment is dehned with endpoints (0.4, 0.4, 0.6)^ and (0.1, 0,0)^. Fifteen 
uniformly spaced points on the line are selected to monitor the accuracy of 
the numerical solutions. 

Finally, we set the true solution of the Robin boundary problem (14.1 p to be 

M(a;) = sinOxsindj/+ 5. (5-11) 


5.2.1 Constant c{x) 

Example 1 c(W) = 1 
In this case, (15. Ill is reduced to 

u{x) = E-{ (5.12) 

which is equivalent to 


roo 

u{x) = E^{ e^^-^°f{Xt)dLt} 


or 


POO 

u{x) = e^^f{Xt)dLt}, 

for a starting point x belonging to the interior of the solution domain. 


(5.13) 

(5.14) 


We will truncate the time interval to [0,T], an approximation to (I5.14|l will 
be 


u(x) = E“ 


Using the fact that 


e‘-f(X,)dLt]. 

(Ax)^ 


dLt ^ {nt - nt-i) 


3e 


(5.15) 

(5.16) 


we can rewrite (15.1511 as 


u(x 


= E^ 




-}• 


3e 


(5.17) 











(a) e = 3Ax, Err = 9.59%, 
NP=1.35e4, Ax=5e-4 

Fig. 4. Cubic domain: number of paths N 
- line segement) 



(b) e = 4Ax, Err =8.84%, 
NP=1.6e4, Ax=5e-4 

2e5 and c(Xt) = 1. (Left - circle; right 


Next identifying the time interval with the length of sample path NP, we have 


u(x 


= 





(5.18) 


where j' denotes each step of the path and j denotes the steps where the path 
hits the boundary. 


At each step along a path we hrst evaluate 




3e 


if Xt. hits the boundary, we then compute f{Xt^){nt. 


np_i) 


(Ax)" 

3e 


followed 


by multiplying it by 3' , which uses the cumulative time of Lt from f = 0 
to tj. Finally, the expectation is done via the average over N sample paths. 


The simulation results of a cubic domain are presented in Fig. 4 and 5. The 
two hgures show the convergency of the approximations as the length of path 
increases from 1.35e4 to 1.43e4 and 1.6e4 to 1.7e4 over the circle and the line 
segment, respectively. Some deviations are seen at the tail in Figure 4(a) and 
among the middle points in Figure 4(b). Meanwwhile, for the spherical and 
ellipsoid domains (Figure 6 and 7), the approximations are better and the 
errors are relatively smaller especially over the line segments, which are below 
3% in Figure 6(b) and Figure 7(b). 


5.2.2 Variable c(x) 

Example 2 c{Xt) = |x|, x is the hrst component of X^ on the boundary. 
Similar to Example 1, we have 











(a) e = 3Ax, Err = 5.50%, 
NP=1.43e4, Ax=5e-4 

Fig. 5. Cubic domain: number of paths N 
- line segement) 



(a) e = 3Ax, Err = 3.96%, 
NP=6e3, Ax=5e-4 



(b) e = 4Ax, Err = 6.49%, 
NP=1.7e4, Ax=5e-4 

= 2e5 and c(Xt) = 1. (Left - circle; right 



(b) e = 3Ax, Err = 1.24%, 
NP=5.5e3, Ax=5e-4 

2e5 and c{Xt) = 1. (Left - circle; 


Fig. 6. Spherical domain: number of paths N = 
right - line segement) 



(a) e = 3Ax, Err = 2.42 %, NP= 
4.5e3, Ax=5e-4 

Fig. 7. Ellipsoid domain: number of 
right - line segement) 



(b) e = 3Ax, Err = 2.44%, 
NP=4.5e3, Ax=5e-4 

N = 2e5 and c(Xt) = 1. (Left - circle; 


u{x) = . (5.19) 

It can be seen that c{Xs)dLs and f{Xt)dLt have the same form, so we can 
handle c{Xs)dLs exactly the same way as f{Xt)dLt. Then, we have 









(a) e = 3Ax, Err = 6.35%, 
NP=1.6e4, Ax=5e-4 

Fig. 8. Cubic domain: number of paths 
right - line segement) 



(b) e = 3Ax, Err = 6.66%, 
NP=1.48e4, Ax=5e-4 

N = 2e5 and c{Xt) = |x|. (Left - circle; 


u{x) = 


,2 

^ eT,Lo 1) 37 /(Xi J (rit. 
i'=o 



(5.20) 


Notice that the term 

gSfc=o(5.21) 

cumulates all the information of c{Xt) with respect to the local time from the 
beginning to the current time. If c{Xt) = |a:|, then 


{ NP . 2 /,2 

^ , (5.22) 

where j' denote each step for the path and j denotes the steps where the path 
hits the boundary. 

Numerical results are shown in Figure 8-10 for a cubic, a spherical and an 
ellipsoid domain, respectively with some adjustment in Ax and NP. Here we 
still have similar results for cube with errors around 6.5%. For the sphere, we 
change Ax to 4e — 4 and there are deviation around the middle in Figure 9(a) 
which may explain the overall error only 6.74% while it performs well over the 
line segment in Figure 9(b) with a smaller error of 3.1%. For the ellipsoid, the 
results are similar as in Example 1 and maintain an error below 4%. 


6 Conclusions and future work 


This paper presents a Monte Carlo simulation method to solve the third 
boundary problems associated with Laplace equations. The idea of simulating 





(a) e = 3Ax, Err = 6.74%, 
NP=6.5e3, Ax=4e-4 

Fig. 9. Spherical domain: number of paths 
right - line segement) 



(b) e = 3Ax, Err = 3.10%, 
NP=6e3, Ax=4e-4 

N = 2e5 and c(Xt) = |x|. (Left - circle; 




(a) e = 3Ax, Err = 3.76 %, NP= (b) e = 3Ax, Err = 1.93%, 

5e3, Ax=4e-4 NP=5e3, Aa:=4e-4 

Fig. 10. Ellipsoid domain: number of paths N = 2e5 and c{Xt) = |x|. (Left - circle; 
right - line segement) 

sample paths of SRBM by the WOS within the strip region shows its efficiency 
and accuracy in estimating local time and evaluating Feynman-Kac formula. 
It should be noted that the cases that g 7 ^ 0 needs further work due to the un¬ 
known exit time out of the sphere at each step. For the Poisson equation, the 
contribution of the source term might be computed as a conditional integral 
|19] . Moreover, the proper truncation of time period is unknown, though it is 
proven that the variance of the approximation increases linearly of T [13] . 


For future work, more flexible domains with local convexity will be considered 
as it relates to the calculations of electrical properties such as the conductivity 
of composite materials where the particle shapes plays an important role ESI. 
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